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ABSTRACT 


Detection of small targets in the presence of noise and sea clutter interference 
presents a formidable task in a radar system design. Conventional radar detection schemes 
such as spectral discrimination and noncoherent integration, have been employed with limited 
success. This thesis investigates an improved target detection scheme suggested by Carlson, 
Evans and Wilson [Ref 1], applicable to search radars, using the Hough transform image 
processing technique. The system concept involves a track-before-detect processing method 
which allows previous data to help in target detection. The technique provides many 
advantages compared to traditional techniques. The improved detectability results from better 
use of old energy and which is equivalent to using a three-dimensional filter matched to the 
target trajectories in addition to the conventional target parameters. The questions answered 
by this thesis concern the effectiveness of the Hough transform in achieving improved radar 
target detection and system detection performance, (i.e., probability of detection and false 
alarm rate as a function of signal to noise ratio) for small, slow targets. System design 
concepts are considered and a full environment simulation including sea clutter and noise is 
implemented to determine the algorithm efficiency and performance in various scenarios. 
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I. INTRODUCTION 


A. BACKGROUND 

Detection of small targets in the presence of noise and sea clutter interference, 
presents a formidable task in a radar system design. A small Radar Cross Section (RCS) 
target (order of 1 m 2 or less) should be detected in the presence of wide dynamic range non- 
Gaussian sea clutter and thermal noise. As we recall from radar theory, radar detectability 
is proportional to the amount of energy returned from the target. Conventional noncoherent 
radar signal processing techniques, compare the range azimuth-elevation cell return against 
a threshold and make a detection decision based upon the results of this comparison. No past 
data is usually saved after the decision process, although some information about the target 
might be contained in the discarded data. Detection improvement in those systems is usually 
achieved by using a high range resolution waveform combined with a high speed scanning 
antenna (60-300 RPM), to perform scan-to-scan integration. The target is extracted from the 
clutter background using the fact that target returns tend to correlate from scan to scan, while 
the clutter returns tend to decorrelate from scan to scan. Spectral discrimination techniques 
using Doppler processing [Ref. 2], are not appropriate for the detection of relatively slow 
targets since the target may be stationary or have zero velocity relative to the local mean 
surface waves. Detection performance improvement is required to enable successful operation 
in the modem environment of reduced RCS targets and restricted Electromagnetic Emission 
Control Policy (EECP), which implies brief periods of radar operation. 

One approach that shows promise uses a high speed digital processor to implement 
a type of processing known as track-before-detect [Ref. 3,4,6,7,23], As presented in greater 
detail in the following chapters, the track-before-detect method uses data taken prior to 
thresholding (or detection) which is processed into feasible tentative tracks on a scan-to-scan, 
or revisit-to-revisit basis. Measures of the track strength, track quality and other track 
parameters are generated for each tentative track, tracks are declared or when appropriate 
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thresholds are exceeded. Detection sensitivity improvement using this method comes from 
avoiding loss of information caused by premature thresholding and discarding data processed 
on prior data scans or revisits. In this way information that takes longer than one scan, or 
revisit, to develop can be exploited. Some approaches to the track-before-detect method were 
offered by Reed [Ref. 3], and Bamiv [Ref. 4], which applied dynamic programming to the 
track-before-detect problem. 

Carlson, Evans and Wilson [Ref. 1] introduced a new method for track-before-detect 
which results in improved target detection and applied the method to an air-traffic control 
problem, namely detecting high speed targets in a clutter free environment. This thesis 
investigates the applicability of this technique for the detection of small low speed targets in 
the presence of noise and sea clutter. The technique combines the raw data from multiple 
scans and translates them into multi-dimensional time range space using the Hough transform. 
The Hough transform which is a curve or a "feature" detector, usually used in image 
processing, translates targets that appear as hits or points in the time range space to curves 
or "features" in the Hough domain. The method uses the Hough transform and the inverse 
Hough transform to extract the detection tracks from the raw radar data. Using this technique 
creates a track defined by the detected curve in the Hough domain and reduces the need to 
revisit a detected target to accomplish the acquisition process. The proposed method 
improves the performance of noncoherent integration from scan-to-scan, from beam-to-beam, 
and from range cell-to-range cell. 

B. OBJECTIVES 

The objective of this thesis, is to investigate track-before-detect signal processing, 
applied to the radar returns of small moving or stationary targets, in an additive noise and high 
sea clutter background. A major goal is to determine the feasibility of using the Hough 
transform to accomplish the track-before-detect function. 

It is expected that extraction of targets using the track-before-detect process will 
enhance target detection since all the energy returned from the target is correlated prior 
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to thresholding. This is equivalent to using a three dimensional filter which is matched 
to the target trajectories, in addition to the conventional target parameters. The questions 
to be answered by this thesis concern the effectiveness of the Hough transform in 
achieving improved radar target detection. The methodology used in this thesis is as 
follows : 

► An efficient algorithm to perform the Hough transform and inverse Hough 
transform was derived; 

► The detection performance (i.e., probability of detection (P d ) as a function of 
Signal to Noise Ratio (SNR) curves, and P d as a function of probability of false 
alarm (P fa ) curves) using the Hough transform process, were determined 
parametricaly; 

► Benchmark performance of an optimal signal processing (using all the target 
energy) was determined; 

► Hough transform processing was compared against the optimum detector in 
Rayleigh interference to determine detection loss; and 

► Hough transform performance as a function of radar and target parameters was 
determined. 

C. THESIS OUTLINE 

The analysis and results associated with the investigation of the above described 
problem are presented in the five subsequent chapters. Chapter II presents the track- 
before-detect concept and the Hough transform technique as an efficient way to 
implement track-before-detect in search radars. In Chapter III a theoretical background 
for radar detection theory of targets in the presence of thermal noise and sea clutter is 
introduced. Chapter IV introduces the algorithm developed in this thesis. The results 
obtained using those algorithms and comprehensive comparison between those results and 
the theoretical results are discussed in Chapter V. Recommendations and conclusions are 
presented in Chapter VI. Finally, the MATLAB computer code for each algorithm is 
presented in the appendix. 
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II. TRACK-BEFORE-DETECT PROCESSING 


A. TRACK-BEFORE-DETECT BASIS 

Track-Before-Detect (TBD) processing [Ref. 21], has a primary objective of 
enhancing target detectability. In contrast to conventional approaches, TBD is an 
algorithm for detecting trajectories, as opposed to detecting single dwell target returns and 
assembling trajectories from those returns. As a result, TBD is applied over multiple 
observation intervals, or dwells, of the radar rather than over the multiple pulse repetition 
intervals that constitute an individual dwell. Consequently, the time between updates in 
the TBD process is typically in the order of seconds. Significant target motion can occur 
between updates making it necessary to associate, over several updates, the resolution 
cells that correspond to the target motions. The name track-before-detect is derived from 
the fact that a sequence of associated resolution cells, comprising a time history of the 
target or track exists before detection is declared. Considering a coherent search radar, 
capable of measuring the range, range rate, and azimuth angle of the target, the track- 
before-detect algorithm consists of the following set of procedures: 

► Integration of the received radar signals over the dwell time of the scanning radar 
beam; 

► Comparing the integrated signals in each resolution cell on each scan to an 
amplitude threshold, which is set lower than normal; 

► Over a number of scans N, forming sets of associated cells, called templates, 
which correspond to possible target motions; 

► In each template count the number of scans for which the amplitude threshold is 
exceeded, calling this number m; 

► Applying m out of M test (i.e., declare target detection if in a template m^M) 
where M is a selected integer less than N; and 

► Continue to form new templates and extend those that have yielded detections by 
repeated application of the rules of association. 
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Figure 1, introduced by Schleher [Ref.5], illustrates the TBD concept implemented 
as a linear velocity filter for short exposure targets. The figure illustrates normalized 
target and sea clutter returns, plotted on a scan to scan basis. The TBD filter is 
implemented as a linear velocity template which is matched to a particular target 
trajectory. One template is required for each realizable target velocity, consistent with the 
radar range resolution. In this example, detection occurs when at least 9 target responses 
fall within the linear velocity template set for a target exposure of 6 seconds (12 scans). 
False alarm occurs when at least nine random sea spikes and noise match the template 
response. To the human eye, a target track is easily visible. However, if only one data set 
is viewed at a time it is impossible to tell with much certainty if a target is present or not. 

In an automatic detection system [Ref. 6], a detection decision can be made by 
thresholding a score that is a function of all the amplitudes along each potential target 
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Figure 1. Track-Before-Detect Process. 
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trajectory, rather then thresholding after each dwell or revisit, and then trying to make 
a decision based on a single data point. A simple scoring function can be defined by 
noncoherent integration along the candidate trajectory, in which case the TBD functions 
as an adaptive noncoherent integrator which follows the target motion. Thus, the TBD can 
be applied in cases where conventional noncoherent integration would be inefficient 
because of target motion. TBD exploits the fact that, because of target motion constrains, 
the position cell of the target return will be correlated from dwell to dwell, whereas the 
positions of the noise spikes will not. 

B. THE HOUGH TRANSFORM 

The Hough transform, (HT) was first introduced by Paul Hough [Ref. 8] in 1961 
as a method of detecting complex patterns in binary image processing. It achieves this by 
determining specific values of parameters which characterize these patterns. Spatially 
extended patterns are transformed, so that they produce spatially compact features in a 
space of possible parameter values. The HT converts a difficult global detection problem 
in image space into a more easily solved local peak detection problem in a parameter 
space. The key idea behind the HT can be illustrated by considering sets of collinear 
points in an image. A set of image points (x,y), which lie on a straight line can, be 
defined by a relation f, such that 

fi(m,c)Ax,y))=y-mx-c= 0, (2.1) 

where m and c are two parameters, the slope and intercept, which characterize the line. 
Equation (2.1) maps each value of the parameter combination ( m,c ) to a set of image 
points. The HT uses the idea that Equation (2.1) can be viewed as a mutual constraint 
between image points and parameter points, and therefore, it can be interpreted as 
defining a one to many mapping from an image point to a set of possible parameter 
values. This corresponds to calculating the parameters of all straight lines, which belong 
to the set that passes through a given image point (x,y). This operation is called back- 
projection of the image point; the definition relation which achieves this is given by 
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g((x,y),(m,c))) =y-xm -c= 0. 


(2.2) 


In the case of a straight line, each image point (x,y) back-projects or defines a straight 
line in the (m,c) parameter space. Figure 2a is a typical point image, and Figure 2b shows 
the parameter lines produced by back-projecting the image points into parameter space 
using Equation (2.2). Points which are collinear in image space all intersect at a common 
point in parameter space, as demonstrated in Figure 2c, where the coordinates of this 
parameter point characterizes the straight line connecting the image points. The HT 
identifies those points of intersection in the parameter space. Determination of the point 
of intersection in the parameter space is a local operation, which is considerably easier 
than detecting extended point patterns in image space. The extension of the method to 
detect parametericaly defined image curves, other than straight lines, is straightforward. 
Image points on a curve, characterized by n parameters a ,..., a , can be defined by an 

i n 

equation in the form of 

/(a iV ..,a n ),(x,y))=0. (2.3) 

By interchanging the roles of the parameters and variables, Equation (2.3) can be used 
to derive the defining relations for the back-projection mapping of image points to 
parameter space 

g((d iV ..,d n ),(x,y))= 0. (2.4) 

This back-projection equation maps out a hypersurface in the n-dimensional parameter 
space. The most probable parameters for image curves are indicated by the intersection 
of several of those hypersurfaces. 
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Duda and Hurt [Ref. 8] suggested that straight lines might be usefully parameterized by 
the length p, and the orientation, 0, of the normal vector to the line from the image origin 
given as 

p =x cos(0)+y sin(0). (2.5) 

This has distinct advantage over the (m,c) parameterization which has a singularity for 
lines with large slopes, i.e., Using the (p,0) parameters means that image points 
map into sinusoidal curves in a two parameter space. 

Once the HT has been accumulated, the pattern of the counts in the accumulator 
is analyzed to estimate the presence and location of the local peaks. The most common 
method is to determine a global threshold. Any accumulator cell with more counts than 
the threshold indicates a possible occurrence of the search for shape. The threshold is 
chosen either using prior knowledge, or it can be automatically selected by analyzing the 
distribution of counts in the array , e.g., a fixed fraction of the maximum count on any 
single accumulator bin. The shortcomings in the HT, as presented by Hunt [Ref. 19], are; 

► The HT suffers from performance degradation due to presence of noise in the 
image plane; 

► The HT tends to favor long lines at the expense of short lines, that is a result of 
the peak detection embedded in the back-projection to the image plane; and 

► The HT implicitly assumes uniform distribution of the noise. 

As a result, image planes that are contaminated with noise, or containing objects with 
almost but not quite parallel edges, will suffer from performance degradation. 

C. USE OF THE HOUGH TRANSFORM IN SEARCH RADAR 

Traditional search radars provide range, azimuth, elevation and Doppler data as 
functions of time, where time is quantized by the scan period, and the coverage sector is 
searched by sequentially looking in all beam directions. Conventional radars usually apply 
noncoherent integration and M-out-of-N decision rule, so that a target is declared if the 
returns for a beam position and a range gate are above a predefined threshold. When 
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considering a moving target that could move from one detection cell to another, the 
integrating efficiency drops; this problem is known as the range-walk problem. The Hough 
transform, as discussed in the previous subsection, is a curve or a "feature" detection method 
which is well suited to locating lines in a plane of noise. When using conventional maritime 
radars, the effective information is usually obtained through time domain analysis,where 
features are mainly due to sea clutter and clutter characteristics which makes the 
discrimination between the target and the clutter a difficult task. Sampling of a radar power 
returns from a target and analyzing them in the time-range plane, leads to a curve which 
represents the target trajectory (where the amplitude of the curve represents the returns 
power level). Detection of this curve reveals all the current information about the detected 
targets. Figure 3 shows an example of a 20 knots radial velocity target in the time-range plane 
for a particular beam and Doppler bin. The slope of the line is determined by the velocity of 
the target. Stationary target, or stationary clutter, would appear as horizontal lines in this 
space. All moving targets would have some finite slope approaching zero, as the target moves 
faster. The time axis runs from the past to the current time and contains the number of scans 
integrated. The target movement, which defines a line in the time-range domain, can be 
defined by the angle 0 of its perpendicular, from the origin and the distance p from the origin 
to the line along the perpendicular, as demonstrated in Figure 3. The Hough transform map 
points in the time-range space into curves in the Hough parameter space (0, p) by 

p =r cos0+/ sin0. (2.6) 

The r and t in this equation are measured from the origin of the time-range space. The 
mapping to the Hough space is done by stepping through 0 from 0° to 180° and calculating 
the corresponding p. Carlson, Evans and Wilson [Ref 1] have shown that by trigonometric 
manipulation. Equation (2.6) can be expressed as 


p =\jr 2 +t 2 sin(0 +atan~). 


(2.7) 
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Figure 3. Time-Range Space Geometry and the Corresponding Hough Parameters. 

Equation (2.7) leads to a sinusoid with amplitude and phase depending on the r-t value of the 
mapped data point. Figure 4 shows the Hough transformation of the data points in Figure 3. 
It can be observed that each p-0 point in the Hough space represents a single straight line in 
the r-t space with the corresponding p-0 value. Each of the curves in the Hough space 
represents the set of all possible lines in the r-t space that cross the corresponding data point. 
The line in the r-t domain that connects all the data points corresponds to the point of 
intersection of all the mapped sinusoid in the Hough space. The time-range space is divided 
into cells equal to the number of range gates times the number of scans being considered. 
Convenient way for mapping the time -range space to the Hough space using simple matrix 
multiplication was offered by Carlson, Evans and Wilson [Ref 1] as follows: 
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Figure 4. Hough Parameter Space Corresponding to Figure 3 Data Points. 

A data matrix RT comprised of radar returns in the time-range space is defined, such that the 
columns contain I radar range returns and the corresponding time value for a particular 
antenna direction 


RT-- 


r \ r 2 r 3 


i 2 



(2.8) 


Now a transformation matrix H is defined as the combination of all the sines and cosines, as 
defined in Equation (2.6) 
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H= 


cos0 sin0 

1 1 

cos0 sin0 

2 2 


cos0 sin0 

N N 


(2.9) 


where 0 takes V values from 0° to 180° using predefined resolution cells. Determination 
of the resolution cells is resolved by the simulation and will be discussed in Chapter IV. 
The multiplication of RT and H gives a N by I matrix R containing the values of p for 
the different data points. The subscripts of p values represent the index numbers of the 
data point transformed and the angle 0 used in the Hough mapping 




R= H*RT= 


( 2 . 10 ) 


P 


i.e 



Each column of the matrix R represents the N p values for each one of the Hough space 
sinusoid. It is clear that the more data points in the image plane we have the bigger R and 
RT will be, resulting in heavier computational load. The transformation matrix is 
precomputed, and its size depends only on the parameter space quantization. 
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III. RADAR DETECTION THEORY 


A. OPTIMUM DETECTION PROCESS 

An optimum receiving process for a return signal of known form with unknown 
phase, which may be used as a reference for comparison of radar detection performance, 
is diagrammed in Figure 5. 


Known Known Time of Desired Pd and 

Spectrum Unknown Phase Peak Signal Pfa 



Figure 5. Optimum Unknown Phase Detection Process. 

All signal parameters, except carrier phase, are known exactly. The receiving filter is 
matched to the signal spectrum. When the signal envelope out of the filter reaches its 
maximum value, the signal to noise average power ratio is equal to the ratio of total 
received signal energy to the noise spectral density. The ratio of signal power at the peak 
of the carrier cycle to rms noise is twice as great and is taken into account in setting the 
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detection threshold. This system will be used to establish .the detection performance of 
a conventional detection process and enables comparison to the Hough based detectors. 

1. Detection Performance of an Optimum Detection Process. 

Consider a sine wave signal present, along with Gaussian noise (with mean 
a and variance a n ), in the input to an IF filter. The resulting output of an envelope 
detector (with amplitude r) is Ricean, i.e., the density function for R is 


/ 


P R ( r ) = — exp| 
a 

n 


\ 


2 2 \ 
r +a 


2a 


' ar^ 


a 


(3.1) 


where I 0 (z) is the modified Bessel function of zero order and argument z. We use the 
transformations 


2 



n 



(3.2) 


where R corresponds' to the signal to noise ratio and y to the normalized output of a 
square law device. Probability of detection is defined as the probability of y exceeding 
a threshold level f, given x and is obtained by the equation 


P\yH\R] = fexp[-(y+R)]I o (2jRy)dy. (3.3) 

In the absence of target, i.e., noise-only case, R=0 and we obtain the probability of false 
alarm as 

oo 

P fa =J exp(-y)dy. (3.4) 

€ 

If instead of a single observation, N observations y n , n=l,2 . . ., N, of the output of a 
square law device with corresponding signal to noise ratio are made and summed, then 

N 

the probability that the sum y exceeds a threshold £ for specified R„ is given by 

t n 

71 = 1 
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N -1 


pj & t ) a p(f>^i *,*,.*„) = 


v JVP' 


exp[-(y-AW)] I 


N -1 


rfy, (3-5) 


where 

A 

R= 

2. Albersheim’s Detection Equation 

Walter Albersheim [Ref. 9] derived a simple formula for signal to noise ratio, 
which is required to achieve a given level of performance for envelope detection of 
nonfluemating signal, in narrow band noise. Despite its simplicity, this empirical equation 
is remarkably accurate. The equation replaces the use of many graphs and cumbersome 
calculations. The Albersheim equation is given as 


N 




,x . 


(3.6) 


[ 4.54 

SNR=- 51og M + 6.2 +- 

10 V ^JM+0.44' 


log io (A +0. 12AB +1 .IB), (dB) (3.7) 


where 


A=ln] 


0.62 




B= Id 


' V 

( 1 -P 


(3.8) 


=ln(P rf )-ln(l-P f ). 


M is the number of independent pulses, P d detection probability and P /a the false alarm 
probability. 

3. Collapsing Loss 

Practical radar systems can seldom preserve their full RF signal through the 
integration and thresholding process. Consider n video samples containing signal plus 
noise integrated along with m extra samples of noise alone. Collapsing loss [Ref. 19] is 
defined as the additional signal to noise ratio required to provide a specific detection 
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performance (P 4 P fa ) for the n+m signal plus noise variates, as compared to that required 
in the case of n signal plus noise variates and is given as 


L = 10 log 


V R , 

1 


(3.9) 


where R, is the required signal plus noise ratio with m extra noise variates and R 2 is the 
signal to noise ratio required with no extra noise variates. The resulting collapsing ratio 
p is defined as 


n+m m 

p = - = 1+—. (3.10) 

n n 

Several empirical methods are available for estimating the collapsing loss. The 
most highly developed method is Batron’s method [Ref. 10], which uses the concept of 
an empirically determined detector loss function. To carry on with the Barton method, it 
is necessary to introduce a quantity called the detectability factor, which is defined as the 
ratio of single pulse energy to noise power per unit bandwidth, that provides stated 
probabilities of detection and false alarm. The notation for the detectability factor is D t (n), 
where i=0, 1, 2, 3, 4 identifies either a steady case (0) or a Swerling type fluctuating 
target (1, 2, 3, 4), while n indicates the number of pulses which are incoherently 
integrated. The detectability factor for a fully coherent detection is D c (n). The empirically 
determined detector loss C n (x) is given by 


C(n) = 


[D o (n)+ 2.3p] 

D ( n) 

o v ' 


(3.11) 


From the definition of the detector loss, the n-th pulse incoherently integrated detectability 
factor is given by 

D (1) C (n) 


DM = 


(3.12) 


n 


Substituting Equation(3.11) into Equation (3.12) and solving for D 0 (n) results in 
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The additional signal to noise ratio required, as compared to coherent integration, can be 
expressed in terms of integration loss given by 



This can be used to find the collapsing loss which is given by 


(3.15) 

A coherent detection curve for D c (1), in terms of P d and P fa is achieved using empirical 
results and are given in [Ref. 10]. 



B. DOUBLE THRESHOLD PROCESS DETECTION PERFORMANCE 

Double threshold, or M-out-of-N detection algorithm, is widely applied in search 
radars. In this process, a fixed group of N pulses is applied to a detector threshold £. A 
detection is declared if the number of threshold crossings within the group exceeds a 
secondary threshold M. For a nonfluctuating target, it was shown that the M-out-of-N 
probability of detection P D , is related to the single pulse probability of detection p d , by 
the binomial relation 
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(3.16) 


k = M' > 

In the absence of a signal, the probability of a false M-out-of-iV detection, due to noise, 
Pfa> is similarly related to the single pulse probability of detection p fa by 

p ^i <3 - i7) 

k = iA > 

The single pulse probability of false alarm is related to the threshold level £ by 

p f = exp(-ri), (3.18) 


where 



n 


(3.19) 


a n is the standard deviation of the noise, and £ the threshold setting. Since it is developed 
in the absence of any signal except for that due to system noise, Equation (3.17) is 
directly valid for all target fluctuation cases (Swerling cases). 

1. Probability of Detection for Nonfluctuating Target 
The nonfluctuating target has been modeled by Marcum [Ref. 11], wherein each 
pulse of a group of N pulses is of constant amplitude. The probability density function 
of such a process is Ricean distributed and may be written in terms of the single pulse 
probability of detection 


oo 



(3.20) 


Equation (3.20) can be recognized as the Marcum Q function. Thus, it follows from 
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Equations (3.16) and (3.20) that for the nonfluctuating case 

pj?) - E f f| (3- 2 i) 

DM k = M \ k 

where Q(R,r\) is the Marcum Q function, and R normalized signal to noise ratio as 
defined by Equation (3.2). 

2. Probability of Detection for Swerling Case I Target 

The slowly fluctuating target has been described by Swerling as one in which 
pulses are perfectly correlated within the integration time but are scan to scan 
independent. For this analysis, the Swerling case I target may be considered to be a series 
of N Marcum model pulses accruing successively but with statistically independente 
amplitudes from scan to scan. If we assume that statistics of this amplitudes variation are 
Rayleigh distributed, we can write 




P m (R> dR, 


(3.22) 


where P DM (s) is the Marcum case detection probability, and G the average power signal 
to noise ratio. Combining Equations (3.21) and (3.22), we obtain R/Jh the detection 
probability for Swerling case I target 


P d £ ft) 0-23) 

/ J ° l 2G I * = jf l k J 

0 

C. DETECTION PERFORMANCE IN THE PRESENCE OF SEA CLUTTER 
1. K Distributed Clutter Model 

For radars in which the dimension of the range resolution cell is much greater than 
the sea swell wavelength with grazing angles greater than about 10°, it is well known 
that the clutter amplitude is Rayleigh distributed [Ref. 12]. This is the consequence of the 
central limit theorem, since the results can be thought of as being the vector sum of two 
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randomly phased components from a large number of independent scatterers. As the radar 
resolution is increased and for a smaller grazing angle, the clutter amplitude is observed 
to develop a longer "tail", and the returns are often described as becoming spiky [Ref. 12, 
13, 14]. It has been found from practical measurements over a wide range of conditions 
that the clutter returns for a high resolution radar can be well modeled by two 
components. The first component is an underlying mean level obeying a chi distribution. 
This mean level has a long temporal decorrelation period, which characterizes for 
example, the mean level variation in clutter spikes, or the periodic variation in amplitude 
as seen when looking up or down the sea swell. The second component is termed the 
"speckle" component. This speckle component has a mean level determined by the first 
component of the clutter model. The K distribution is derived by averaging the speckle 
components over all possible values of the mean level 

oo 

p(r) =J p(r/u)p(u)du, (3.24) 

o 


where p(r) is the overall Probability Density Function ( PDF) of the clutter returns, p(u) 
is the PDF of the clutter mean level, and p(r/u) is the PDF of the speckle component 
.The speckle component is Rayleigh distributed so that 


ixr 

P(r/u) = — exp| 
2 u 


2 \ 


■nr 


\ 


4 u 


r^O, 


(3.25) 


where u obeys a chi distribution and p(u) is given by 

P(u) = K 2v-1 exp(-fc V), 0, (3.26) 

T(v) 

where T is the gamma function, v is a shape parameter, b 2 = v/Q and Q=E[u 2 ] is the 
average power of the clutter. Substituting Equations (3.25) and (3.26) into Equation (3.26) 
yields the K-distribution 


= 4 T^y rV ^-i (2cr) ’ 


(3.27) 
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where K v (x) is the (v-l)-th order modified Bessel function of the third kind, and c=b/^J2 

is the distribution scale parameter. For high resolution sea clutter, values of v are 
generally observed in the region 0.1<v<°o, where v~ 0.1 represents a very spiky and v=°° 
represents thermal noise. An empirical derived shape factor for an x band radar is offered 
by Watts [Ref. 12] and given by 


log(v)=log 


2 5 

- log(cj) )+— 

J 8 O 


(3.28) 


where <}> g is the grazing angle in degrees, l is the cross range resolution in meters, i[j is 
the aspect dependency, which takes value of for up or down swell directions, — for 
across swell and 0 for intermediate directions, or no swell, and k describes the 
polarization effect and takes the value 1 for vertical polarization and — for horizontal 
polarization. 

2. Detection Performance in K Distributed Clutter 

The analysis for the detection performance in K-distributed clutter [Ref. 15, 16] 
assumes that the detection processing consists of a linear bandpass filter and a square law 
detector, followed by a pulse to pulse integrator. This receiver is not optimal for the Re¬ 
distributed clutter, but it is widely used in practice, since the detection is easy to 
implement. The output of a square law detector is sampled in turn at time intervales t=iT 
(i=0,land then integrated. A test statistic V is compared with a threshold rj to 
decide whether or not a target exists. The output of the receiver can be represented by 

r(t) = S(t)+C(t)+n(t), (3.29) 

where S (t), C(t ) and n(t) represent signal, K-distributed clutter and receiver noise, 
respectively. They can be expressed as 
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S(t) = a(t) cos(o t)-b{t) sin(o) t), 

c c 

C(t ) = Jt^(t)cos(oW) -y(f)sin(co^f), (3.30) 

n(t ) = x (f)cos(w ?)-v (?)sin(co t), 

n c n c 

where a(t), b(t), x£t), x n (t) and y n (t) represent the inphase and the quadrature component 
of the signal, clutter, and noise, respectively. They are mutually independent random 
processes. Equation (3.29) can be rewritten as 


r(t) =[a{t) +x(f)]cos(<of) ~[b(t) +y(f)] sm(co^), (3.31) 

where 


x(t) = x (f)+x (t), 

c n 

y(t) =c(t)+y (f). 


(3.32) 


are the inphase and quadrature of the clutter plus noise. The average power of x{t) or y{t) 
is given by 


E[x 2 (f)] = E[y\t)] = f(u 2 +p 2 )p(u)du = — + P 2 = Q + P 2 , 

J n *2 n n 

0 u 


(3.33) 


where P n 2 is the average power of white noise component, and Q is the average power 
of the clutter. The Average Clutter to Noise Ratio (CNR) is given by 

v Q 

Vi 1 = J 1 ' (334) 

n n 


A = 


The test statistic can be defined as 


l w-i 

V - — E Vx+af^y+bA 

2P 2 t - o 11 11 


(3.35) 


where and b, (i=Q,l,...,N-l) are the samples of x(t), y(t), a( t) and b(t), respectively. 
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For the convenience of the analysis, u in Equation (3.25) will be regarded as a random 
variable. Based on this assumption, the conditional PDFs of the inphase and quadrature 
components of the clutter plus noise can be expressed as 


P(x/y)=- 


N 


2tx(m 2 +P 2 ) 

n 


p(y/u)=- 


\ 


2tt(h 2 +P 2 ) 

n 



(3.36) 


The Moment Generation Function ( MGF), defined as the Laplace transform of the PDF 
of the statistic V in the presence of nonfluctuating target, can be obtained by [Ref. 10] 


L y (zJu,R) = 
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(3.37) 


In the case of a fading target the R in (3.37) is a random variable. For a chi square family 
of fluctuating targets the PDF of R is given by 


p(R) = 


1 




m\ r) 


k ( 

B k -1 

R exp 


kR 

R' 


(3.38) 


where R is the average of R over target fluctuations and k > 0 is a fluctuating parameter 
representing the degree of signal fading. The smaller k corresponds to a deeper signal 
fading, k can be related to the Swerling’s fluemating model, where k=l, N, 2, 2N and k=°° 
corresponds to Swerling cases 1, 2, 3, 4 and constant target, respectively. The MGF of 
test statistic V of the chi square family of fluctuating target can be written as 
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L y (zju) = J~.L y (zJu,R) p(R)dR. 
o 


Solving the integral leads to 
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(3.39) 


(3.40) 


In the absence of target, i.e., R= 0, the MGF becomes 

1 


L o (zJu) = 



(3.41) 


The conditional false alarm probability can be obtained using the residue technique to 
solve the inverse Laplace transform 


p M = -£res 
fa w 


exp(riz) 

L(z/u) -, z, 


kO 


(3.42) 


where res [•] denotes residue, r| is a decision threshold and z k0 (ko= 1,2,...) are the poles 
of L 0 {z/u) lying in the left half plane. Solving Equation (3.42) results in 


T| n expl 
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(3.43) 
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The false alarm probability is then obtained by integrating all possible values of u 
obeying chi-distribution in Equation (3.26) 
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In the case of K>N, 

u 2 NR 

J 2+ T 

n 

' nr) n ~" 

k) 

The detection probabilities are obtained by averaging Equations (3.48) and (3.49) over all 
possible values of u that is 

oo 

P d =fp d (u)p(u)du. (3.50) 

o 

where p(u) is given by Equation (3.26). 
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IV. ALGORITHM DESCRIPTION 


A. SYSTEM CONFIGURATION 

This chapter discusses the system implementation and the design considerations 
of the various parameters used in the simulation. The system block diagram is presented 
in Figure 6. 



Figure 6. System Configuration. 

















The algorithm is comprised of three main blocks. The first block is the radar data 
simulation, responsible for generating the time range date map. This block accepts the 
following inputs: 

► Radar system parameters, such as Pulse Width (PW), antenna scan rate, peak 
power and system losses; 

► Environmental conditions such as sea state, radar grazing angle and sea swells 
characteristics. Those parameters determine the clutter behavior, as discussed in 
Section B; and 

► Target characteristics, which define the target radar cross section, Swerling case 
type, location in space and maneuvering information. 

The output of this function is a random variable generated as the coherent sum of the 
three components of the radar returns, namely target signal, noise and sea clutter. The 
output from the first block serves as the input to the Hough based detector function. This 
block performs the task of the Hough transform, inverse Hough transform and system 
detection performance evaluation. Conventional detection schemes performance evaluation 
is done at the third part of the system and establishes a reference baseline for detection 
performance comparison. 


B. RADAR DATA SIMULATION 

To simulate a maritime radar data map, we need to generate the three main 
components of the radar returns, i.e., target signal, noise and sea clutter. We generate 
signal plus interference random variable given by 



2 2 
( vs-rn ) +rn , 

i 9 


(4.1) 


where vs is the steady signal voltage, rn { and rn q are the inphase and quadrature 
components of the interference respectively. Two types of targets are considered in this 
work; the non fluctuating target and the Swerling type I target. The nonfluctuating target 
is modeled as a constant signal, whose amplitude is determined by the radar equation 
given as 
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where P t is the radar peak power, G the antenna gain, a the target cross section, R the 
detection range and L the overall system losses. To simulate a Swerling I target, we use 
a model [Ref. 18] which describes the fluctuation in target amplitude caused by the 
changes in target aspect angle, rotation or vibration of target scattering sources, or 
changes in radar wavelength. The Swerling case 1 model assumes a Rayleigh target 
amplitude fluctuation, which results in an exponential cross section distribution given by 


- 1 , 
w(o,o)-—exp 
a 


a 


v 


(4.3) 


where a is the target mean radar cross section. 

The clutter addressed in this work is sea clutter and considered to be K distributed. 
A method of simulating K-distributed clutter introduced by Schleher [Ref. 17], suggests 
that the K-distributed random variable is given by 


r = u^J-2 ln(x), (4-4) 

where x is a uniformly distributed variable (0-1), u is root gamma (chi) variable defined 
as 



(4.5) 


where b is a scale parameter, defined by Equation (3.34), and t is a gamma random 
variable whose PDF is given by 


Pit ) = 


r(v) ’ 


tz 0. 


(4-6) 


Equation (4.6) implies that K-distributed clutter simulation requires the generation of a 
gamma random variable of order v. This is done by first generating a gamma deviate 
of integer order k and then a second gamma deviate of order v-k, such that, a= v-k< 1. 


The sum of the two gamma variates provides a gamma variable of order v. A gamma 
variable of integer order k can be achieved by adding up k exponentially distributed 
waiting times, i.e., the logarithms of uniform deviates given by 

ti=-M x** 2 * *x k ), (4.7) 

where x. are uniformly distributed random variables (0-1) and tl k =0. The second Gamma 
random variable of order v-k is generated using the rejection method, which is a general 
technique for generating random variables whose PDF is known and computable. Figure 
7 demonstrates the geometrical argument on which the rejection method is based. 



to 


Figure 7. Rejection Method for Generating Random Deviate t with a Known PDF 
P(t). 
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We consider a known PDF, p(t) denoted as the shadowed region in Figure 7. The first 

step in the rejection method consists of selecting an arbitrary invertible function f(t), 

defined as the comparison function, which lies everywhere above the original probability 

function. Next consider a point in two dimensions that is uniform in the area under the 

comparison function (points a and b in Figure 7). If the test point lies outside the original 

probability function, the point will be rejected; if it lies inside the area, it will be 

accepted. The accepted points are uniform in the accepted area, so their t values have the 

desired distribution. To choose a random point in two dimensions, a test random t is 

generated from the comparison function, using a uniform deviate ranging from zero to A, 

where A is the area under the comparison function. A second uniform deviate (0-1) is 

p(t) 

generated and tested against the ratio q= -. If the deviate is less than q, the variable 

t is accepted; otherwise it is rejected. An efficient comparison function for the gamma 
random variable is 


M 


m 

0<tel 

-t 


c 

tel 

l r(a) 



(4.8) 


Random variables t can be generated from /(f) using the transformation method 


Kx)=F~\t), (4.9) 

where F ! (t) is the inverse function of the indefinite integral of f(t). Applying Equation 
(4.8) to (4.9) and solving results in 
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aT(a)x a 
al\a)(l -ln(x)) 


0<t<l 
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The rejection criteria is then obtained by computing 




0<t<l 

tel 


which leads to 


(4.10) 


(4.11) 
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At this point, receiver noise is introduced into the simulation. The noise amplitude is 
assumed to be Rayleigh distributed; and therefore, the noise power amplitude is 
exponentially distributed (we assume a square law envelope detector). Since the noise 
envelope is independent of the clutter, its variance adds to that of the clutter. The 
resulting random variable simulating k-distributed clutter with receiver noise is given by 




(-2 ln(x)). 


(4.12) 


where P„ is the noise power. Substituting Equation (4.12) into Equation (4.1) and (4.4), 
yields the target plus interference random variable given as 


rns= 




(vs-m cos(2tcjc j )) 2 +(rn sin(2nx i )) 2 , 


where x 1 is uniform distributed random variable. 


(4.13) 


C. HOUGH TRANSFORM BASED DETECTOR 

The Hough transform algorithm transforms data from the time-range space to the 
Hough space, using the formulation of Chapter III. Figure 8 presents the configuration of 
the time-range space and the Hough space parameters used for the simulations. The time 
range space is divided into 128 range gates and a time history composed of 100 scans, 
with scan time defined by the radar parameters. This leads to time-range space dimensions 
of 128 by 100 for a total of 12800 data cells. 
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Figure 8. Time-Range Space and Hough Space Configuration. 


Figure 9a is a time-range mesh plot example of an approaching constant target with a radial 
velocity of 20 knots and heading of 45°, relative to the line of sight. Figure 9b is the same 
target with additive Rayleigh noise. The raw data is compared against a primary threshold 
whose level is determined by the requirement that the initial probability of false alarm, which 
is applied to the data at this point, must allow most of the target returns through with some 
noise spikes. The main objective of the primary threshold is to reduce the computation load, 
which is linearly proportional to the number of time-range cells crossing the primary 
threshold. The primary threshold level, rq, is determined by the desired value of the initial 
probability of false alarm (per pulse) in the time-range domain, P fa , given (for the case of 


35 




additive Rayleigh noise) by 




(4.14) 


A gray scale plot of the time-range space for a target signal, with additive noise 
before and after applying a primary threshold corresponding to Pfal of 10' 2 , is shown in 
Figures 10a and 10b, respectively, where only the region of interest is shown, i.e., the 
range cells in which the target exists. 




Time (sec) Range 


(b) 

Figure 9. Mesh Plot in Image Space for a 20 Knot Target. 

(a) Without Noise. 

(b) With Additive Rayleigh Noise (SNR=2dB). 
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Figure 10. Time-Range Space for 20 Knot Target with Additive Rayleigh 

Noise (SNR=5 dB). 

(a) Before Primary Threshold. 

(b) After Primary Threshold. 

The time-range data is now transformed to the Hough parameter domain using the 
procedure described in Chapter III. The signal power associated with the sinusoidal curve 
segments in each cell is now summed noncoherently to give a total cell power. Figure 1 la 
shows the Hough parameters mesh plot corresponding to the data presented in Figure 9b. The 
major peak matches the target trajectory while the smaller peaks arc related to the noise; each 
peak produces a straight line in the time-range space. Applying a secondary threshold in 
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1. Binary Integration Hough Transform 

Implementing the Hough transform when more than one target is present in the 
observed area might result in dynamic range problem [Ref. 1], If one target is larger than the 
other, we get masking of the small target by the large target in the Hough space. This results 
from the Hough transform integration approach using the (which performs a noncoherent 
summation of the signal power), whereby, a large target peak in the Hough domain will be 
much higher than the small target peak. Setting the threshold low enough to detect the small 
target will increase the number of false alarms, while increasing the threshold results in losing 
the small target. As a result choosing a threshold that allows detection of small targets, while 
avoiding detection of false targets, might be a difficult task. Similar problems occur if we 
detect targets that have crossing trajectories; as a result the peak of one target will mask the 
other target peak. One possible way, discussed by Carlson, Evans and Wilson [Ref 1], to 
avoid this problem is to use binary integration in the Hough space. In this method, all the data 
points that cross the primary threshold in the time-range space are assumed to have a value 
of one and mapped into the Hough space accordingly. This eliminates the difference between 
the targets due to their size and results in a target size in the Hough space relative to the 
targets trajectory length. Numerous conventional detection schemes which employ a fixed 
threshold, or use a M-out-of-N secondary threshold might be used to detect the target peaks. 
Using the binary integration method has some disadvantages, since noise or clutter peaks that 
cross the primary threshold level might cause a false alarm. 

2. Optimum Parameter Selection 

Several parameters can affect the performance of the Hough based detector. The first 
one is granularity of the parameter space. The optimum granularity is the one that achieves 
the desired P fa without choosing too fine a granularity that splits the target between more 
than one cell in the Hough domain. This issue was investigated by applying different 
granularity to the Hough transform and observing the effect on the system performance. 
Other parameters that should be addressed are the primary and secondary threshold value 
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selection. We can choose high primary threshold values, which result in fewer data points 
crossing the threshold and less points transformed into the Hough domain. In this case, we 
have to use low values for the secondary values, as the target trajectory is comprised of less 
points. The other alternative is to use low primary threshold values, which will cause most of 
the time-range points to be transformed into the Hough domain and determine high secondary 
threshold values to avoid false alarms. Finding the optimum values for detection is done by 
using different values of threshold in the simulations and evaluating their detection 
performance. Those results are presented in Chapter V. 

D. HOUGH TRANSFORM DETECTION PERFORMANCE 

Efficient evaluation of the Hough transform based detectors, and a measure of 
comparison to the more conventional detectors, is achieved by determining the detection 
statistics (i.e., probability of detection as a function of the input SNR for a given P fa ). This 
is done for specific radar parameters and for predefined target, noise and sea clutter signal 
levels. The detection statistics will enable us to evaluate the target detection range in a given 
scenario. 

1. Probability of False Alarm 

The first step in evaluating the detector performance consists of determining the 
probability of false alarm as a function of secondary threshold. False alarm in the Hough 
detector is defined as finding a line, or a trajectory, in the Hough domain where no target 
exists in the time-range space. To find the probability of false alarm, we use the concept of 
"accessible Hough space" introduced by Carlson, Evans and Wilson [ Ref. 1]. This quantity 
is found by putting a one in all the range time cells of the image space, and then Hough 
transforming them. The resulting mesh plot is shown in Figure 12 for the Hough parameters 
defined in Section B . It can be observed that some cells in the Hough domain have a zero 
value which results from the fact that not all the cells in the Hough domain can be 
backprojected into the time-range domain , in other words those cells match non-existing 
trajectories. “ The accessible Hough space” is parameter dependent and has to be calculated 
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according to the granularity of the time-range and the Hough dimensions. To determine the 
Pfc as a function of secondary threshold, a time-range map with Rayleigh noise was generated 
using the simulation procedure described in chapter V. The primary threshold in the time- 
range domain is set to a predefined value (several values were tested), and the Hough 
transform is applied to this data. The number of secondary threshold crossings, divided by the 
total number of cells in the accessible Hough space, yields the probability of false alarm. The 
secondary threshold is then varied to determine the probability of false alarm, as a function 
of the secondary threshold. 


Rho 


0 


Theta 



Figure 12. Plough Transform of All Ones in Image Space, demonstrating the 
"Accessible Hough Domain". 
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2. Probability of Detection 

Detection for the Hough based detector is defined as finding a predefined target 
trajectory in the time-range space. To evaluate the detection performance of the system, we 
first define a scenario and then test it under those conditions. In the scenario definition, we 
simulate a target trajectory in the time-range space, and then combine the simulated Rayleigh 
noise and sea clutterwith the target according to the desired SNR. Transformation to 
parameter space then gives the corresponding accumulator cell vaiues in Hough space. 
Probability of detection is calculated by counting the number of detections achieved in the 
appropriate Hough cell divided by the overall number of simulations. As was proposed by 
Carlson, Evans and Wilson [Ref. 1], an n by n grid about the expected Hough cell location 
is considered, where n is an integer defined in the simulation, and detection is declared as any 
secondary threshold crossing within this grid. The proper value for n is determined in the 
simulation. The procedure is repeated for different secondary threshold values, corresponding 
to different values of P fa , which will give the system probability of detection as a function of 
the SNR for a given P fa . 




V. RESULTS 


A. TEST CASES DESCRIPTION 

The simulated radar parameters, specified in Table 1, are typical for land based 
high resolution radars. These types of radar are mainly used for small target detection in 
the presence of sea clutter and, as such, are potential candidates for the Hough detector 
implementation. The algorithm was applied to different scenarios generated in accordance 
with the radar parameters, chosen to enable comprehensive evaluation of the system 
performance and efficiency. The parameter space for this simulation is taken as 256 by 
180 for the p and 0 dimensions. Those parameters, as shown in subsection E, are close 
to the optimum choice for the granularity of the parameter space. 


Table 1. Simulated Radar Specifications. 


Center Frequency 

9 GHZ 

Pulse Repetition Frequency (PRF) 

1000 

Pulse Width (PW) 

10 nsec 

Peak Power 

10 MW 

Antenna Gain 

35 dB 

Antenna Speed Table 1. 

240 RPM 

System Losses 

6 dB 

Grazing Angle 

2 

Polarization 

Horizontal 


Two main test scenarios are considered in this work. The first one consists of a single 
target, nonfluctuating or Swerling 1 type, embedded in Rayleigh noise or K-distributed 
sea clutter. Several target speeds, heading directions and SNR’s are considered to 
determine the system sensitivity to those parameters. The second scenario consists of 
multiple targets with additive noise and K-distributed sea clutter, with variation in the 
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targets relative size, speed and headings. This scenario is intended to evaluate the 
performance of the algorithm in which multiple targets present in the time-range space. 

B. SINGLE TARGET SCENARIO 

1. Additive Rayleigh Noise Case 

A simple but important test case is a single constant velocity target with additive 
Rayleigh noise. The target is tested for different speeds and heading directions, where 
those parameters define the target trajectory in the time-range space and affect the system 
performance. The different test cases for this scenario are summarized in Table 2, where 
R,, denotes the target range, V 0 is the target velocity, ao the target acceleration and a 
denotes the targets heading direction relative to the radar site. The SNR is assumed to be 
constant during the integration time, due to the low target speed and the resulting small 
crossing distance. The primary threshold in the time-range space was chosen to give a 
time-range P fa of — and the secondary threshold value, £, which is taken as a fraction 
of the maximum value of the accumulator cell in the parameter space. 


Table 2. Single Target Test Cases 


Test Case I 

Stationary Target 

Nonfluctuating 

Ro=20 Nm, V o =0 

Test Case II 

Constant Velocity 

Target 

Swerling I type 

Ro=20 Nm, V 0 =25 Knots, 

a =90° 

Test Case III 

Accelerating Target 

Swerling I type 

Ro=20 Nm, V 0 =5 Knots, 
a =45° a^O.l ft/sec 2 


Figures 13a-14a are the time-range space maps after primary threshold for test case I, 
with SNR of 0 and 5 dB, respectively, and £ matched to the peak value. Figures 13b-14b 
present the true target trajectory and the Hough detector time-range back-projection, while 
Figures 13c-14c show the corresponding Hough parameter space. The longest possible 
trajectory for this particular time-space map is 162, occurring at the diagonal line, while 
the maximum target trajectory length, which corresponds to this test case, is 100. We can 
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observe the false detection declared in the case of SNR=0 dB. These occur because the 
detector prefers the longer lines created by the noise. In the case of 5 dB SNR, the peak 
in the parameter space is higher, and the target is easily detected. The effect of choosing 
an unmatched threshold is shown in Figures 15-16, where the optimal threshold is reduced 
by 1 dB and applied to the test case I target, with SNR’s of 0 and 5 dB, respectively. The 
resulting false tracks demonstrates the sensitivity of the detector to the secondary 
threshold level. 

Test case n has a diagonal trajectory in the time-range space, which corresponds 
to the longest line possible in the time-range map. Figures 17-18 show the time-range 
space map and the corresponding parameter space for SNR’s of 0 and 5 dB, respectively, 
with matched £. Again we can observe the clear peak in the Hough domain corresponding 
to the target trajectory. The effect of unmatched £ is tested again and is presented in 
Figures 19-20 for secondary threshold setting of 0.5 and 1 dB lower than the optimal one. 
Test case m consists of a target changing its speed by a radial acceleration over the data 
history. In order to detect this type of target, we need to map the data into a 3-D Hough 
space. A 2D detector, like the one implemented in this thesis, can be used, and the results 
of applying the algorithm to test case III are shown in Figures 21-22 for SNR of 2 and 
7 dB, respectively. A lower threshold level is used, and we can see that several lines are 
mapped back into the image space. Each of those lines represent a section along the 
actual curved path. The true trajectory is tangent to those lines along the inside edge. 
Although some uncertainty results from using a 2D system to detect a 3D space, it is not 
extreme, and it is much easier to implement. 

As discussed in Chapter IV, the Binary-Integration (BI) Hough detector has some 
advantages in the multiple target scenarios and might be used as a preferred detection 
scheme. To evaluate the effectiveness of this scheme, the BI Hough detector was applied 
to test cases 1-0, and the results are shown in Figures 23-27 for SNR’s of 0 and 5 dB and 
matched £. The detection performance reflected from those plots seems to be lower than 
those obtained using the Hough detector itself, although for SNR’s of 5 dB the original 
target trajectory is being detected in all the test cases. Quantitative results for the 
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performance degradation using this method, as compared with the results using the basic 
algorithm, are presented in Section D. 

2. Additive K Distributed Clutter Case 

A more realistic scenario is the case of one target in the presence of sea clutter. 
The clutter was simulated in accordance with the procedure described in Chapter IV, 
where the clutter characteristics are simulated using Equation (3.28), in accordance with 
the radar parameters and assuming up/down swell direction. We tested the algorithm 
against test cases I-II for different interference levels where interference power is defined 
as I=P +P =P (1 +CNR), P c and P n are the noise and clutter power, and CNR defines the 

c n n 

clutter to noise ratio. Figures 29-30 present the results obtained for test case I with 
additive K-distributed clutter for signal to interference ratio (SIR) of 0 and 5 dB, 
respectively, and matched £. Similar results to the additive noise case are obtained for 
thi s case, although the clutter has spikier nature than the noise. We can observe that the 
target trajectory is detected in the case of SIR of 5 dB, while no detection is declared in 
the case of 0 dB. 

The results for test case 2 are presented in Figures 31 and 32 for SIR’S of 0 and 
5 dB, respectively, with matched £. The effect of unmatched threshold is demonstrated 
in Figures 33-36 for test cases I-II and the two tested SIR’s. The overall performance of 
the detector in the presence of K-distributed clutter is close to the performance in the 
noise-only case, so that the effect of the sea spikes, which affect the performance of 
conventional detectors, has been eliminated. 
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Figure 15. Test Case I, SNR=0 dB, £=18.5 dB. 

(a) Time Range Space after Primary Threshold. 

(b) True and Detected Target Trajectory in Image Space. 

(c) Parameter Space Mesh Plot. 











Figure 16. Test Case I, SNR=5 dB, ^ = 18.5 dB. 

(a) Time Range Space after Primary Threshold. 

(b) True and Detected Target Trajectory in Image Space. 

(c) Parameter Space Mesh Plot. 
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Figure 18. Test Case II, SNR=5 dB, £=21 dB. 

(a) Time Range Space after Primary Threshold. 

(b) True and Detected Target Trajectory in Image Space. 

(c) Parameter Space Mesh Plot. 
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Figure 21. Test Case III, SNR=2 dB, £=16 dB. 

(a) Time Range Space after Primary Threshold. 

(b) True and Detected Target Trajectory in Image Space. 

(c) Parameter Space Mesh Plot. 
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Figure 22. Test Case III, SNR=7 dB, £=20 dB. 

(a) Time Range Space after Primary Threshold. 

(b) True and Detected Target Trajectory in Image Space. 

(c) Parameter Space Mesh Plot. 
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Figure 24. Test Case I, SNR=5 dB, 5=19 dB. 

(a) Time Range Space after Primary Threshold. 

(b) True and Detected Target Trajectory in Image Space. 

(c) Parameter Space Mesh Plot (with Binary Integration). 
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Figure 26. Test Case II, SNR=5 dB, £=18 dB. 

(a) Time Range Space after Primary Threshold. 

(b) True and Detected Target Trajectory in Image Space. 

(c) Parameter Space Mesh Plot (with Binary Integration). 
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Figure 27. Test Case IE, SNR=5 dB, £=16 dB. 

(a) Time Range Space after Primary Threshold. 

(b) True and Detected Target Trajectory in Image Space. 

(c) Parameter Space Mesh Plot (with Binary Integration). 
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Figure 28. Test Case I, with Additive K-distributed Clutter, SIR=0 dB, £=22 dB. 

(a) Time Range Space after Primary Threshold. 

(b) True and Detected Target Trajectory in Image Space. 

(c) Parameter Space Mesh Plot. 
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Figure 30. 

(a) 

(b) 

(c) 


Test Case II with Additive K-distributed Clutter, SIR=0 dB, £ =23 dB. 
Time Range Space after Primary Threshold. 

True and Detected Target Trajectory in Image Space. 

Parameter Space Mesh Plot. 
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Figure 31. Test Case II, with Additive K-distributed Clutter, SIR=5 dB, t,= 22 dB. 

(a) Time Range Space after Primary Threshold. 

(b) True and Detected Target Trajectory in Image Space. 

(c) Parameter Space Mesh Plot. 
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Figure 32. Test Case III, with Additive K-distributed Clutter, SIR=5 dB, £=22 dB. 

(a) Time Range Space after Primary Threshold. 

(b) True and Detected Target Trajectory in Image Space. 

(c) Parameter Space Mesh Plot. 
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Figure 33. Test Case I with Additive K-distributed Clutter, SIR=0 dB, £ = 16 dB. 

(a) Time Range Space after Primary Threshold. 

(b) True and Detected Target Trajectory in Image Space. 

(c) Parameter Space Mesh Plot (with Binary Integration). 
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Figure 34. 

(a) 

(b) 

(c) 


Test Case I, with Additive K-distributed Clutter, SIR=5 dB, £=17 dB. 
Time Range Space after Primary Threshold. 

True and Detected Target Trajectory in Image Space. 

Parameter Space Mesh Plot (with Binary Integration). 
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Figure 35. Test Case II with Additive K-distributed Clutter, SIR=0 dB, £=18 dB. 

(a) Time Range Space after Primary Threshold. 

(b) True and Detected Target Trajectory in Image Space. 

(c) Parameter Space Mesh Plot (with Binary Integration). 
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Figure 36. Test Case II, with Additive K-distributed Clutter, SIR=5 dB, £=19 dB. 

(a) Time Range Space after Primary Threshold. 

(b) True and Detected Target Trajectory in Image Space. 

(c) Parameter Space Mesh Plot (with Binary Integration). 
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C. TWO TARGETS TEST SCENARIO 

As discussed in Chapter IV, one of the potential problems using the Hough based 
detector occurs when integrating time-range maps containing more then one target. Two 
multiple target scenarios are considered to determine the performance in those cases. The 
first one, test case IV, defines two targets with crossing trajectories and RCS ratio of 1:4 
between the targets. Figure 37 presents the time-space and the corresponding parameter 
space for test case IV with SNR of 7 dB. Although the target trajectories are similar in 
length, we can observe the size of the different peaks in the parameter domain resulting 
from the targets relative size. Due to the high SNR in this test case and with an apriori 
knowledge about the existence of two targets in the integration area, we can set the 
secondary threshold, so that the two targets are detected. When we consider the same 
scenario with low SNR, the small target is masked by the noise, and it is impossible to 
detect it using a single threshold. This phenomena is demonstrated in Figure 38 for a SNR 
of 2 dB. Applying the BI Hough detector to the test case IV targets are presented in 
Figures 40-41 for SNR’s of 7 dB and 2 dB, respectively. This results in the emphasis of 
the small target in the parameter space, and an equal peak size in the parameter space for 
the two targets is observed. Applying a simple single threshold in this case enables 
detection of the two targets. 

The second case of interest is the case of a small target masking by a very large 
target. This is defined for the case of high SNR in test case V, which defines two 
crossing targets with RCS ratio of 1:10 between the targets. Figure 41 shows the time- 
range domain and the resulting Hough parameter for test case V with SNR of 15 dB. In 
spite of the fact that the two target trajectories are very close in length, the large target 
peak masks the small target, since it is much higher in power; and as a result the small 
target can not be detected. The BI Hough detector was applied to this case, and the results 
are shown in Figure 42. The level of the two targets peaks in the Hough domain are 
almost equal now, which enables detection of the two targets using a single threshold. 
However, the BI Hough transform has some performance degradation for low SNR’s, 
which is discussed in Section C. 
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Figure 37. Test Case IV, SNR=7 dB, £ = 16 dB. 

(a) Time Range Space after Primary Threshold. 

(b) True and Detected Target Trajectory in Image Space. 

(c) Parameter Space Mesh Plot. 
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Figure 39. Test Case IV, SNR=7 dB, £ = 16 dB. 

(a) Time Range Space after Primary Threshold. 

(b) True and Detected Target Trajectory in Image Space. 

(c) Parameter Space Mesh Plot (with Binary Integration). 
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Figure 40. 

(a) 

(b) 

(c) 


Test Case IY, SNR=2 dB, £=12 dB. 

Time Range Space after Primary Threshold. 

True and Detected Target Trajectory in Image Space. 
Parameter Space Mesh Plot (with Binary Integration) . 
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Figure 41. Test Case V, SNR=15 dB, £=20 dB. 

(a) Time Range Space after Primary Threshold. 

(b) True and Detected Target Trajectory in Image Space. 

(c) Parameter Space Mesh Plot. 
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D. DETECTION PERFORMANCE 


1. Additive Rayleigh Noise Case 

The first step in performance evaluation is the determination of the probability of 
false alarm as a function of the secondary threshold. We simulate Rayleigh noise in 
image space, with image and parameter space dimensions, as defined in Section A. The 
number of accessible Hough cells, which is determined by applying all ones in the time- 
range space, was found to be 8192 for this particular grid choice. The primary threshold 
was then set to 1/200, which is chosen to reduce the computational load. To obtain 
accurate results for the 5000 simulations were used for each secondary threshold point, 
and the results were averaged. Figure 43 presents the results for the probability of false 
alarm as a function of secondary threshold. 

To find the probability of detection, a single target was included in the simulation. 
The test target is modeled as a nonfluctuating target, approaching the radar site with 
different speeds, corresponding to a diagonal line in the time-range space. To obtain 
accurate results, the algorithm was tested with different SNR’s and 500 repetitions per 
SNR point. Figures 44-46 show the resulting detection curves for target speeds of 0, 5, 
and 25 knots, respectively. The different curves in those plots represent different 
secondary threshold settings corresponding to different P fa ’s derived from Figure 43. 
Figures 48 -50 show the corresponding detection curves for the BI Hough algorithm with 
P fa , as a function of the secondary threshold derived from Figure 37. Observing the 
detection curves of the BI Hough algorithm, we find about 1 dB performance degradation 
relative to the Hough detector. The Hough detector performance is compared against the 
optimum detection process and binary integration detection process performance. Figure 
51 presents the P d as a function of SNR for a given P fa for a stationary target with a P fa 
of 3-10 4 . Figure 52 presents the performance comparison for target speeds of 5, 15 and 
25 knots. Figures 53-54 show the performance comparison for the same scenario with P fa 
of 1-10" 7 . The best detection performance for the conventional schemes are achieved, as 
expected, when the target is stationary. When the target is moving, collapsing loss 
degrades the detection performance of those detectors. The Hough detector performance 
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improves as the target speed increases, resulting in longer target line in the time-range 
domain. Figure 55 shows significant results obtained by estimating the processing 
gain/loss of the detector, which is compared to the optimal detector for different speeds. 
This figure shows the range-walk problem exists when applying traditional detection 
techniques. The Hough detector does not suffer from this problem, and its performance 
improves as the target speed increases. 

2 . Additive K-Distributed Clutter Case 

Detection performance in the case of additive K-distributed clutter is evaluated in 
a similar way to the noise only case. The probability of false alarm as a function of the 
secondary threshold is shown in Figure 56, where the time-range and Hough space 
parameters are as defined in Subsection 1. Comparing the results obtained in this case to 
those presented in Figure 43 for the noise-only case, we observe that the P fa , or the same 
secondary threshold setting, is higher, due to the spiky nature of the clutter. 

The probability of detection was determined using the test targets defined in 
Subsection 1. The number of repetitions in this case was limited to 100 per SNR point, 
due to the large computational load in the K-distributed clutter simulation. Figures 57-59 
present the resulting detection curves for a scenario consisting of target embedded in k- 
distributed clutter plus noise approaching the radar site in velocities of 0, 5, and 25 knots, 
respectively. The different curves in those plots represent different secondary thresholds. 
Comparison of the Hough detector to traditional radars is done by simulating a square law 
detector and evaluating its performance for the same scenarios. Figure 61 presents the 
performance comparison for a nonfluctuating stationary target, a 25 knot nonfluctuating 
target and a P fa of 8-10" 5 . We observe that the Hough transform performance improves as 
the target speed increases. The spikier nature of the clutter, which force us to use a higher 
threshold to achieve predefined P fa level, results in causes degradation in the performance 
comparing to the noise only case. Figure 60 presents the processing gain/loss for this 
case, as compared to an envelope detector for different target speeds. We observe that the 
range-walk problem is solved in this case, as in the noise-only case although the 
processing gain is lower in higher velocities (2 dB). 
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Probability of False Alarm 



Figure 43. Hough Detector P fa as a Function of Secondary Threshold. 
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Probability of Detection 



Figure 44. Hough Detector P d as a Function of SNR for a Nonfluctuating Stationary 
Target with Additive Rayleigh Noise Scenario. Different Curves Denote 
Different secondary Threshold Settings. 
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Probability of False Alarm 



Figure 45. Hough Detector P d as a Function of SNR for a Nonfluctuating 5 Knot 
Target with Additive Rayleigh Noise Scenario. Different Curves 
Denote Different secondary Threshold Settings. 
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Probability of Detection 



Figure 46. Hough Detector P d as a Function of SNR for a Nonfluctuating 25 Knot 
Target with Additive Rayleigh Noise Scenario..Different Curves 
Denote Different secondary Threshold Settings. 
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Probability of Detection 



Figure 48. Binary Integration Hough Detector P d as a Function of SNR for a 

Nonfluctuating Stationary Target with Additive Rayleigh Noise Scenario. 
Different Curves Denote Different secondary Threshold settings. 
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Figure 49. Binary integration Hough Detector P d as a Function of SNR for a 

Nonfluctuating 15 knot Target with Additive Rayleigh Noise Scenario. 
Different Curves Denote Different secondary Threshold Settings. 
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SNR (dB) 


Figure 50. Binary integration Hough Detector P d as a Function of SNR for a 

Nonfluctuating 25 knot Target with Additive Rayleigh Noise Scenario. 
Different Curves Denote Different secondary Threshold Settings. 
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Probability of Detection 



Figure 51. P d as a Function of SNR for a Nonfluctuating stationary Target 
with Additive Rayleigh Noise scenario @ P fa =3-10 -4 . Different 
curves denote different detection schemes. 
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Probability of Detection Probability of Detection 



(a) 



(b) 

Figure 52. P d as a Function of SNR for a Nonfluctuating Target with Additive 

Rayleigh Noise scenario @ P fa =3-10 -4 , Different curves denote different 
detection schemes. 

(a) Target Speed 5 Knots. 

(b) Target Speed 25 Knots. 
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Probability of Detection 



Figure 53. P d as a Function of SNR for a Nonfluctuating stationary Target 
with Additive Rayleigh Noise scenario @ P fa -1-10' 7 . Different 
curves denote different detection schemes. 
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Probability of Detection Probability of Detection 



(a) 



(b) 

Figure 54. P d as a Function of SNR for a Nonfluctuating Target with Additive 

Rayleigh Noise scenario @ P fa =l-10' 7 . Different curves denote different 
detection schemes. 

(a) Target Speed 5 Knots. 

(b) Target Speed 25 Knots. 
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Required SNR (dB) Required SNR (dB) 




Figure 55. Required SNR to Achieve P d =0.9 as a Function of Target Speed in 
Knots, for Optimum Detection Process and the Hough Detector. 

(a) P fa = 3-10 4 . 

(b) P fa = HO' 7 . 
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Probability ot Detection 



Figure 57. Hough Detector P d as a Function of SNR for a Nonfluctuating Stationary 
Target with Additive K-Distributed Clutter Scenario. Different Curves 
Denote Different secondary Threshold Settings. 
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Probability of Detection 



SIR (dB) 


Figure 58. Hough Detector P d as a Function of SNR for a Nonflactuating 15 knot 
Target with Additive K-Distributed Clutter Scenario. Different Curves 
Denote Different secondary Threshold Settings. 
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Figure 59. Hough Detector P d as a Function of SNR for a Nonflactuating 25 knot 
Target with Additive K-Distributed Clutter Scenario. Different Curves 
Denote Different secondary Threshold Settings. 
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Probability of Detection Probability of Detection 



P d as a Function of SNR for a Nonfluctuating Target with Additive 
K-Distributed Clutter Scenario @ P fa =8- 10 s . Different curves 
denote different detection schemes. 

(a) Stationary Target. 

(b) Target Speed 25 Knots. 




















Required SNR (dB) 



Figure 61. Required SNR to Achieve P d =0.9 @ P fa = 8-10' 5 as a Function of Target 
Speed in Knots, for Square Law Detection Process and the Hough 
Detector, (for the Additive K-Distributed Clutter Scenario). 
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E. PARAMETERS SELECTION 

1. Parameter Space 

To determine the optimum granularity in the Hough space, we consider a 
nonfluctuating 25 Knot constant velocity target using the additive Rayleigh noise scenario. 
Using the procedure described in Section D we evaluate the detection performance, i.e., 
P d as a function of required SNR, for the different granularity in the parameter space. 
Figure 62a shows the required SNR to achieve a P d of 0.9, with secondary threshold of 
12 dB (corresponding to P fa of 110' 7 ). The 0 granularity was fixed at A 0=1° and Ap 
values were varied between 0.25° to 2.75°. Figure 62b presents the required SNR to 
achieve P d = 0.9 for the same scenario with fixed Ap and A0 ranging from 0.25° to 
2.75°. The results obtained for the different granularity cases follow the general expected 
behavior. Selecting a very fine granularity results in splitting the target into several 
parameter cells, while a coarse granularity causes the noise to add up in fewer cells, 
forming false peaks in the parameter space, which dictates higher values of £ to achieve 
the same P d . The optimum choice is close to the granularity hat we used in the simulation, 
i.e., A0=1° and Ap=l. 

2. Primary and Secondary Thresholds 

The second factor that affects the Hough detector performance is the primary and 
secondary threshold settings. As discussed in Chapter IV, the P fa is determined by both 
the primary and secondary thresholds. To determine the optimal values for each, we 
evaluate the P d as a function of required SNR for different values of primary threshold 
and apply secondary threshold s that yields constant P fa . Figure 63 shows the results for 
the required SNR to achieve P d of 0.9 and a P fa of 1-10' 7 , as a function of the £ for the 
scenario of nonfluctuating 25 knot constant velocity target with additive Rayleigh noise. 
The primary threshold is varied from —to We observe that the optimum choice for 
the primary threshold is about —. The primary threshold used to evaluate the detection 
performance is which is a compromise between optimum detection performance and 
the required computing power when applying a low primary threshold. 
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1.5 

Rho Granularity 



1 1.5 

Theta Granularity 


Figure 62. Required SNR a Function of Parameter Space Granularity 
£=12 dB. 

(a) A 0=1, p varies. 

(b) A p=l, 0 varies. 


Pd=0.9, 
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Required SNR (dB) 



Figure 63. Required SNR as a Function of Primary Threshold @ Pd=0.9, £=12 dB. 
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VI. CONCLUSIONS 


In this thesis, the track-before-detect processing method, using the Hough transform 
for target detection and tracking by a high resolution search radar, was examined. An efficient 
algorithm for Hough detector implementation offered by Carlson, Evans and Wilson [Ref. 1], 
was used, and a comprehensive simulation was built to test the algorithm for several 
scenarios. The system performance was evaluated and compared against traditional detection 
schemes. The detection capability of the detector for small targets in the presence of noise and 
sea clutter, was proven, as well as, the tracking feature obtained by the system as part of the 
detection process, which eliminates the need to revisit the target for tracking information. 

The algorithm presents an improved detection performance, compared to conventional 
detection schemes, when the scenario is composed of moving targets embedded in Rayleigh 
noise or K-distributed sea clutter. This provides a solution to the range-walk problem, which 
occurs in traditional detection schemes where targets with inter-cell motion are not 
completely integrated. The problem is solved using this method, since the algorithm detects 
trajectories rather than single resolution cells. A scenario with a 25 knot constant velocity 
target in additive noise, shows a processing gain of about 3 dB, relative to the optimal 
stationary detector for this particular parameter selection. When considering a stationary 
target with additive Rayleigh noise, the detection performance is significantly degraded 
compared to the optimal detector, or to the binary integration detection process. The required 
SNR to achieve P d =0.9 with a P a of 10' 7 is about 5 dB for this case. This results mainly from 
the limited integration time that the detector uses and the resulting relatively short trajectories 
in the time-range domain. The detector employed was designed to detect straight trajectories 
in the time-range domain. However, the algorithm was able to find the tangent to the curves 
trajectory of an accelerating target, which indicates the feasibility of using this algorithm in 
more complex scenarios than a constant speed target. Evaluating the detector performance 
in the presence of K-distributed sea clutter, which is the main interest of maritime search 
radars, shows promising results and the processing gain for a 25 knots constant velocity 
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target was found to be 1.5 dB. 

The main drawback of the method is its sensitivity to parameter selection, both in the 
time-range space and in the Hough space. The radar parameters must be well matched to the 
target's kinematics to enable satisfactory results. Using unmatched radar range resolution, or 
scan rates, causes the target to stay in the same resolution cell and reduces the detector 
performances dramatically, as no significant line in the time-range domain exists. The same 
performance degradation results from improper selection of parameter space granularity. 
Using course granularity leads to noise adding up to significant lines while fine granularity 
results in the target splitting into several parameter cells. Optimization evaluation for the 
Hough space granularity reveals that the optimum values for A0 and Ap should be 
determined in accordance with the time-range and Hough space parameter selection. Another 
trade off that should be considered when applying the algorithm is the one between primary 
and secondary threshold settings. The results obtained in this work show that an optimal 
choice between the two thresholds exists. However, choosing a low primary threshold level 
implies higher computational load; therefore, the threshold setting should be selected in 
accordance with the available computing power. This might dictate using higher than 
optimum primary threshold level to enable real time application of the algorithm. A possible 
detector problem occurs for multiple target scenarios when one target masks the other 
targets, or when small targets are masked by noise. One optional solution uses the binary 
integration Hough transform, which allows multiple target resolving; however, the detection 
performance in this case is degraded by about 1 dB. Other alternatives to solve this problem 
involve multiple threshold algorithms and were not investigated in this thesis. 
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APPENDIX A . MATLAB SOURCE CODE 


The algorithms described in chapter IV were implemented using MATLAB software. The 
different functions are listed below. 


A. HOUGH TRANSFORM 

function [RT,Q,angnuxn]=ht(bi,type.dat) 

% Filename : ht.m 
% Title: Hough Transform 
% Date of last Revision: 10 Mar 95 
% Comments: 

% This m file computes the Hough transform using the classical H- Transform. 
% Input variables: 

% bi- flag which indicates if we perform HT (0) or binary integration HT (1) 

% type- Defines the environment simulation noise(O) or noise plus clutter (1) 

% dat-File name contains the radar, target and space parameters 

% Output variables: 

% RT- Matrix containing the Hough transform values of the radar echoes. 

% Q- Matrix contains the power of the accumulative Hough transform 
% angnum- Number of required angle bins in the Hough transform. 

% Associated Matlab Functions: 

% rtmap.m (generates time-range map) 

% iht.m 

% Associated Matlab files 
% rdata.m (radar parameters) 

% Initialization 
eval([dat]); 

[trmap,pwr,mgres,pmap,echo ) T,dt ] =rtmap(dat,type); 
angnum= 181/angres; 

N=linspace(0,180,angnum); 
tl=0:dt:T; 
t=0:length(tl)-l; 
if bi==l 

pwr=ones(size(pwr)); 

end 

% Hough transform 
H2=[sin(N’*pi/180)]; 

H1=[cos(N’*pi/180)]; 
dimen=size(trhit) 
sumlen=0; 
for p=l:dimen(2) 
len(p)=max(size(find(trhit(:,p)))); 
for q =l:len(p) 

R( 1 :length(N),sumlen+q)=H 1 *trhit(q,p)+H2*t 1 (p); 
end 


105 




sumlen=sumlen+len(p); 

end 

if sum(len)~=0 

offset=abs(min(R(:))); 

end 

RT=(round((R+offset.*ones(size(R)))./rhores))+l; 
if ( max(RT(:)))>rhonum 
RT=RT.*rhonum ,/max(RT(:)); 
end 

dimh=size(RT); 

Q=zeros(angnum,max(RT(:)); 
for l=l:dimh(2) 
for k=l:dimh(l) 

Q(k,RT(k,l))=Q(k,RT(k,l))+pwr(l); 

end 

end 

B. INVERSE HOUGH TRANSFORM 

function [rho,tet,mg]=iht(Q,T,dt) 

% Filename : iht.m 
% Title: inverse Hough Transform 
% Date of last Revision: 10 Mar 95 
% Comments: 

% This m file computes the inverse Hough transform 
% Input variables: 

% Q- Matrix contains the power of the accumulative Hough transform 

% T- Integration time in time-range domain 

% dt - Time interval between processed dwells 

% Output variables: 

% rho- The rho value where the peak in the Hough domain was detected 
% tet- The theta value where the peak in the Hough domain was detected 
% mg- Range locations of the detected trajectory 
% Associated Matlab Functions: 

% ht.m 

% trmap.m 

% Associated Matlab files 
% rdata.m 

t=0:dt:T; 

snr2=max(Q(:))/0.913; 

[tet,rho]=find(Q>snr2); 
ntr=length(rho); 
for test=l:50 
if nrt>3 

snr2=snr2* 1.0233; 
end 
end 

for k=l:length(rho); 
rho 1 (z)=(rho(k)-1 )*rhores-offset; 
r0(k)= rhol(k)/cos(tet(k)*pi/180); 
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a(k)=-tan((90-tet(k))*pi/180); 
mg(k,:)=t./a(k)+rO(k); 
end 

C. RADAR DATA PARAMETERS 


% Filename : rdata.m 
% Title: radar parameters 
% Date of last Revision: 10 Mar 1995 
% comments: 

% m file that contains the simulated radar and target parameters 
% Input variables: 

% None 

% Output variables: 

% Radar and target parameters as defined in the program 
% Associated Matlab Functions: 

% ht.m 

% iht.m 

% trmap.m 

% Associated Matlab files 
% None 


% Radar and environment parameters 
rpm=240; % antenna speed 


sigml=10; 

sigm2= 10; 

pt=le6; 

gt=le4; 

lambda=0.1; 

scan=100; 

pw=0.02; 

thetag=0.5; 

rc= 100; 

sl=1/3; 

kl=1.7; 

% Targets parameters 

tetal=0; 

teta2=0; 

r01=10; 

r02= 10.005; 

v01=20; 

v02=25; 

snrdb=20 

num=T; 


target 1 cross section (m A 2) 
target2 cross section (m A 2) 
radar peak power (w) 
radar antenna gain 
radar wavelength 
# of antenna scans used 
radar pulse width 
Grazing angle 
cross range resolution 
swell direction 
polarization factor 


tetal=0; % targetl heading ( relative to the radar) 

teta2=0; % target2 heading ( relative to the radar) 

r01=10; % targetl initial range 

r02= 10.005; % target2 initial range 

v01=20; % targetl speed 

v02=25; % target2 speed 

snrdb=20 % SNR in dB 

num=l; % number of targets to be simulated 

% Time-range and Hough space parameters 

mgnum= 128; % number of range cells to be integration in the time-range space 

rhores=l; % defines the granularity for rho in the Hough domain 

angres=l; % defines the granularity for theta in the Hough domain 

rhonum=256; % defines number of the cells in the transform 

pfa= 1/100; % defines the requires primary P fa in the time-range domain 

bi=l; % flag indicates if we perform HT (0) or binary integration HT (1) 
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type= 1; % Defines environment simulation noise(0) or noise plus clutter (1) 

D. TIME-RANGE MAP GENERATION 


function [trmap,pwr,mgres,pmap,echo,T,dt ] =rtmap(dat,type) 

% File name: rtmap.m 
% Title: Time-Range Map Generation 
% Date of last revision: 10 Mar 1995 
% Comments : 

% This m file creates simulated time range map composed of target noise and clutter 
% Input variables: 

% dat- File name contains the radar and target parameters 
% type- Defines the environment simulation noise(O) or noise plus clutter (1) 

% Output variables: 

% T - Total time length in which the radar echoes are processed 
% dt - Time interval between processed dwells 
% pwr -Vector containing the power of the received radar echoes. 

% trmap - Matrix containing the hits crossing primary threshold (2 by length(t)) 

% rangres - The radar range resolution (in ft) 

% pmap- Vector contains the time-range hits power 

% Associated matlab functions 
% ht.m 

% iht.m 

% Associates matlab functions 
% rdata.m 

eval([dat]); 
dt=60/rpm; 

T=scan*dt; 

t=0:dt:T; 

rangres= pw* 150*3.281; 
snr=(10r(snrdb/10); 

[r 1 ,r2,P 1,P2]=tgt(t,dat); 

P=(Pl+P2)/2; 

means=mcan(P); 

sigmn= sqrt(means/(2*snr»; 

thers=-sigmn*log(pfa); 

pmap=zeros(scan,mgnum); 

for i=l:scan+l 

pmap(i,r(i)+1)= pmap(i,r l(i)+1)+(P l(i)); 
pmap(i,rl(i)+1)= pmap(i,r2(i)+ l)+(P2(i)); 
end 

vmap=sqrt(pmap); 
rtot=zeros(length(t),mgbin); 
if type==0 

for i=l:min(size((vmap))) 

[noise]=rayn(mgbin); 

noisei= (sigmn.*noise(:,l))’; 

noiseq= (sigmn.*noise(:,2))’; 

echo=sqrt((vmap(i,:)-noiseq). A 2+noisei. A 2); 

loc=find((echo(i,:))>thers); 
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pwr(i,:)=zeros(size(temp3)); 
pwr(i,loc)=echo(i,loc); 
rtot(i,:)=(loc,zeros( 1 ,mgbin-length(loc))]; 
end 

elseif type==l 
pwr=[]; 

for i=l:min(size((vmap))) 
inrf=kcltr(means,rngbin,dat); 
echo=sqrt((M(i,:)-inrf(:,2)’). A 2+(mrf(:,l)’). A 2); 
loc=find((echo(i,:))>thers); 
pwr=[pwr echo(loc)]; 
rtot(i,:)=(loc,zeros(l,mgbin-length(loc))]; 
end 
end 

trhit=[rtot’]; 

E. TARGET SIMULATION 

function [rl ,r2,P 1 ,P2]=tgt(T,dt,dat) 

% File nam e: tgt.m 
% Title: target simulation 
% Date of last revision: 10 Mar 1995 
% Comments : 

% This m file creates simulated targets in accordance with the radar equation 

% Input variables: 

% dat - File name containes the radar and targets parameters 

% T - Total time length in which the radar echoes are processed 

% dt - Time interval between processed dwells 
% Output variables: 

% rl,2- Target range 

% PI,2- Target received power 

% Associated matlab functions 
% trmap.m 

% Associated matlab files 
% rdata.m 

eval[dat]) 
t=0:dt:T; 

rl=r01*6076-(v01*6076/3600)*t*cos(tetal); 
r2=r02*6076-(v02*6076/3600)*t*cos(teta2); 
r 1=floor((rl ,/mgres)))+0.5); 
r2=floor((r2./mgres)))+0.5); 

Pl=((pt*gt A 2*labbda A 2*sigm)./((4*pi) A 3.*r 1/3.218). A 4»; 

P2=((pt*gt A 2*labbda A 2*sigm)./((4 *pi) A 3. *r2/3.218). A 4)); 
if num-=l 

Pi=[);n=[I; 

end 

F. RAYLEIGH NOISE SIMULATION 

function [ray]= rayn(N) 
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% File name: rayn.m 
% Title: Rayleigh noise Simulation 
% Date of last revision: Mar 10 1995 
% Comments : 

% This m file simulates Rayleigh noise 
% Input variables: 

% N- Number of noise points to be simulated 
% Output variables: 

% ray- Matrix containing the I and Q components of the noise (2xN) 

% Associated Matlab functions: 

% trmap.m 

% Associated Matlab Files: 

% None 

n=l:N; 

mdl =randn(size(n)); 

md2=randn(size(n)); 

rnd=sqrt(rndl. A 2+md2. A 2); 

mdang=2*pi.*rand(size(n)); 

ray=[md’.*sin(mdang)’,md’.*cos(mdang)T; 

G. K-DISTRIBUTION CLUTTER SIMULATION 

function [cltr]=noisel(sig,N,dat) 

% Filename : kcltr.m 
% Title: k distributed clutter simulation 
% Date of last Revision: 10 Mar. 1995 
% comments 

% This m file simulates distributed clutter 

% Input variables: 

% sig-Signal power 

% Output variables: 

% cltr-Matrix containing the I and Q components of the interference (2xN) 
% Associated Matlab functions 
% trmap.m 

% Associated matlab files 
% rdata.m 

eval(dat); 

sf= 10 A (0.67*logl0(thetag)+0.625*logl0(pw* le3*rc/30)+sl-kl); 

pc=sig*(10 A (cnr/10))/(l+10 A (cnr/10)); 

pn=sig A 2/(l+(10 A (cnr/10))); 

f=floor(sf); a=sf-f; 1=0; 

p=exp(l)/(a+exp(l)); 

c=(a+exp(l))/(exp(l)*gamma(a+1)); 

cnt=floor(2*c*s); 

b=sqrt(2*sf/pc); 

ran=rand(N,f); 

for i=l:N 

x(i)=-log(prod(ran(i,:)); 

g(i)=x(i); 

end 
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if a~=0 
for i=l:cnt 
u=rand; v=rand; 
if u<p 

y=v A (l/a); q=exp(-y); 
else 

y=l-log(v); q=y A (a-l); 
end 

if rand<=q 
1-1+1; g(l)=x(l)+y; 
end 

if 1==N 
break 
end 
end 
i=l:N; 

chi(i)=(sqrt(g(i)))./b; 

inrf(i)=(sqrt((chi(i).) A 2+0.5*pn)).*sqrt(-2*log(rand(size(i)))); 
mdang=2*pi.*rand(size( inrf)); 
cltr=[( inrfl.*sin(mdang))’,( inrf.*cos(mdang))T; 

H. HOUGH DETECTOR PROBABILITY OF DETECTION 

function [pd ]=prdet (pl,p2,sl,s2 ,simlen,bi,type,dat) 

% Filename :prdet.m 

% Title: probability of detection of the Hough detector 
% Date of last Revision: 10 Mar 1995 
% Comments: 

% This m file computes probability of detection (P^ of the Hough detector as a 

% function of the secondary threshold , for different primary threshold settings 

% (in steps of 0.25 dB) 

% Input variables: 

% pl,p2- Primary threshold boundaries (in dB) 

% sl,s2- Secondary threshold boundaries (in dB) 

% bi- Flag which indicates if we perform HT (0) or binary integration HT (1) 

% simlen- Number of repetitions in each SNR point ( primary and secondary ) 

% type- Defines the environment simulation noise(0) or noise plus clutter (1) 

% pfa- Defines the requires primary P fs in the time-range domain 
% Output variables: 

% pd- Matrix containing the pd 

% val - Matrix containing the primary and secondary threshold used for pd 

% Associated Matlab functions: 

% ht.m 

% iht.m 

% Associated Matlab files: 

% none 

m=0,n=0; 

for snr2db =sl:0.25:s2 
m=m+l; 

snr2= 10 A (snr2db/10); 
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for snrdb= pl:0.25:p2 
n=n+l; 
hitl=0 

for count=l:simlen 
[RT,Q,angnum]=ht(bi, type.dat); 

[rho,tet,mg]=iht(Q,t); 
for i=l:length(tet) 

if tet(i)<=reftet(tgt)+l & tet(i)>=reftet(tgt)-l 

if rho(i)-offset<=reffho(tgt) +1& rho(i)-offset>=refrho(tgt)-l 
hit(i)=l; 
end 
end 
end 

if sxnn (hit)>0 
hitl=hitl+l; 
hit=0; 
end 
end 

pd(m,n)=hit 1 ,/simlen; 
hitl=0; 
end 
end 

val=[sl:s2;pl:p2]; 

I. HOUGH DETECTOR PROBABILITY OF FALSE ALARM 

function [pfa ]=prfa (sl,s2,simlen,bi,type,dat) 

% Filename : prfa.m 

% Title: probability of false alarm of the Hough detector 
% Date of last Revision: 10 Mar 1995 
% Comments: 

% This m file computes probability of detection (P fa ) of the Hough detector as a 

% function of the secondary threshold (in steps of 0.25 dB) 

% Input variables: 

% sl,s2-Secondary threshold boundaries (in dB) 

% simlen-Number of repetitions in each SNR point ( primary and secondary ) 

% dat-File name contains the radar and target parameters 
% Output variables: 

% pfa-Matrix containing the probability of false alarm 
% val-Vector containing the secondary threshold used for the p ta 

% Associated Matlab functions: 

% ht.m 

% iht.m 

% Associated Matlab files: 

% none 

m=0 ; 

for snr2db =sl:0.25:s2 
m=m+l; 

snr2= 10 A (snr2db/10); 
for count=l:simlen 
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[RT,Q,angnum]=ht(bi, type.dat); 

[rho,tet,mg]=iht(Q,t); 

hit(count)=length(rho); 

end 

pfa(m)=sum(hit)/simlen; 

end 

val=sl:s2 ; 

J. ALBERSHEIM METHOD 

function [snr,pd]=alber(M,N,DC,pfa) 

% File name: alber.m 
% Title: albersheim metod 
% Date of last revision: Mar 10 1995 
% Comments : 

% This m files implement the albersheim equation-required snr to 
% obtain a given pd and pfa including collapsing loss 
% input varaibles: 

% M-Overall number of Signal cells 

% N-Number of overall cells 

% pfa-Required probability of false alarm 

% DC- Detectablity factor for a coherent process ( vector size 1x99) 
% Output variables: 

% snr-Required snr to obtain a given p d and p fa 
% pd- pd corresponds to the snr 
% Associated Matlab Functions: 

% None 

% Associated Matlab Files: 

% None 


pd=0.01:0.01:0.99; 

A=log(0.62/pfa); 

B=log(pd./(l-pd)); 

snrl=-5*logl0(N)+(6.2+4.54/(sqrt(N+0.44))).*logl0(A+0.12.*A.*B+1.7.*B); 

rho=N/M; 

dc=10. A (DC./10); 

1= 10*logl0((l+sqrt(l+9.2*rho*M./dc))./(l+sqrt(l+9.2*M./dc»); 
snr=snrl+l; 


K. BINARY INTEGRATION THRESHOLD 

function [q]=binther(N,M,PFA) 

% File name:binther.m 
% Title: binary integration threshold 
% Date of last revision: Mar 10 1995 
% Comments : 

% This m file computes the required threshold level needed to obtain overall Pfa 
% input varaibles: 

% M- Overall number of Signal cells 
% N-Number of overall cells 
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% PFA- Overall P FA required for the binary integration process 
% Output variables: 

% pfa-Single pulse probability of false alarm 
% Associated Matlab Functions: 

% binsnr.m 

% Associated Matlab Files: 

% None 

steps=-0.1; 
test= 100; 
x= le-9; 
k=M:N; 

delt= le-6*PFA; 
for kk=M:N 

noverk(kk-M+ l)=prod(l :N)/(prod(l :N-kk)*(gamma(kk+1))); 
end 

for loop= 1:500 
p=x/k; 

q=(l-x)/(N-k); 
r=sum(noverk.*p *q); 
if abs(r-PFA)<=delt 
break 
else 

if sign(test)~=sign(r-PFA) 
steps=-0.5*steps; 
end 
end 

x=x+steps; 
test= r-PFA; 
end 

if loop==500 

input (’no convergence in Pfa’.’s 1 ) 
break 
end 

q =-log(x); 

L. BINARY INTEGRATION PROBABILITY OF DETECTION 

function[pd] =binpd(N,M,PD) 

% Title: binary integration probability of detection 
% Date of last revision: Mar 10 1995 
% Comments : 

% This m file computes the requires single pulse p d needed to obtain overall P D 
% input varaibles: 

% M-Overall number of Signal cells 

% N-Number of overall cells 

% PD-Overall P D required for the binary integration process 
% Output variables: 

% pfa-Single pulse probability of detection 

% Associated Matlab functions: 

% binsnr.m 


114 





% Associated Matlab Files: 

% None 
for l=l:length(PD) 
test= 100; 
x=le-4; 
steps=-0.1; 
k=M:N; 
delt= le-6; 
for i=M:N 

noverk(kk-M +l)=prod(i+l:N)/ prod(l:(N-i)); 
end 

for loop= 1:500 
p=x. A k; 

q=(l-x). A (N-k); 
r=sum(noverk. *p. *q); 
if abs(test-PD(l) )<=delt 
break 
else 

if sign(test)~=sign(r-PD(pp)) 
steps=-0.5*steps; 
end 

x=x+steps; 

end 

clear p q 
test= r-PD(pp) ; 
end 

if loop==500 

input (’no convergence in PD’/s 1 ) 
break 
end 

pd(r)=x; 

clear r k i noverk 
end 

M. BINARY INTEGRATION DETECTION PERFORMANCE 

function[SNR]=binpdnr(N,PD,PFA) 

% File name: binsnr.m 

% Title: binary integrationprobability of detection 
% Date of last revision: Mar 10 1995 
% Comments : 

% This m files computes the required snr to obtain a given pd and pfa 
% using binary integration processing 

% input varaibles: 

% M- Overall number of Signal cells 
% N-Number of overall cells 

% PD- Overall P D required for the binary integration process 
% Output variables: 

% snr- Required SNR to achieve to defined P D 
% Associated Matlab functions: 
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% binther.m 
% binpd.m 

% Associated Matlab Files: 

% None 

M=0.4*N; 
test= 100; 
delt=0.0001; 
steps=-.l; 
snr=le-2; 

[q]=binther(N,M,PFA); 

[pd]=binpd(N,M,PD); 
for count=l:length(pd) 
for loop= 1:1000 
h(l)=exp(-snr); 
f(l)=q *exp(-q ); 
d(l)=h(l); 
for k=2:50 
f(k)=q *f(k-l)/k; 
d(k)=snr*d(k- l)/(k-1); 
h(k)=h(k-l)+d(k); 
end 

Q(count)=l-sum(f.*h); 
if abs(Q(count)-pd(count))<=delt 
break 
else 

if sign(test)~»sign(Q(count)-pd(count)) 
steps=-0.5*steps; 
end 

snr=snr+steps; 

end 

if loop==500 
inputfno convergence’,’s’) 
break 
end 

test= Q(count)-pd(count); 
end 

SNR(count)=snr; 
clear d h f Q 
end 


N. GRAPHICS UTILITY 

% Filename : graplot.m 
% Title: plotting utility 
% Date of last Revision: 10 Mar 1995 
% Comments: 

% This m file plots the time-range space and the Hough space parameters 
% Several plotting options are available to choose from. 

% input varaibles: 
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% None 

% Output variables: 

% None 

% Associated Matlab functions: 

% ht.m 

% iht.m 

% alber.m 

% binsnr.m 

% Associated Matlab files: 

% rdata.m 

close all 

k=menu (’Select Required Plot”,Raw Data Mesh Plot’,Raw data Time-Power Plot’,Range -Time Plot’,’Hough 
X-form Mesh Plot’, Trough X-form Rho-Theta Plot*); 
if k==l 
figure(l) 
yl=tl-T; 
xl=0:mgnum; 

subplot(21 l),surf(x 1 ’,y 1,10*log 10(echo)) 
ylabel(Time (sec) 1 ) 
xlabeKRange 1 ) 
zlabel(Rower (db)’) 
elseif k==2 

stairs(tl-T,10*logl0(pmap)), grid 

axis([-T 0 (min(10*logl0(pmap)))-10 (max(10*logl0(pmap)))+10 ]) 

xlabel(’time’) 

ylabel(’power (db)*) 

elseif k==3 
hold on 
plot(Rl,tl-T,Y) 
plot(R2,tl-T,’r’) 
for i= 1 :length(rho 1) 
plot(mg(i,:),tl-T ,’b’) 
end 

xlabeKRangeO.ylabelCTimeO.grid 
title([Target echos range res =’, num2str(rangres),’ft’l) 
elseif k==4 

xl=0:angnum-l; 

yl=min(R(:)):rhores:max(R(:)); 

surf(xl’,yl,Q’),grid 

xlabel(’Theta’) 

ylabel(Rho’) 

zlabel(’Power’) 

title( Hough space parameters’) 
elseif k==5 

N=linspace(0,180,angnum); 
plot (N,R,’r’),grid 
xlabel(’theta’) 
ylabelCRho 7 ) 

title( Hough space parameters’) 
elseif k==6 
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plot(10*logl0(SNR),PD,’b’),grid 


hold on 

plot(10*logl0(snr),pd,’r T ) 

xlabelfSNR (dB)’),ylabel(Pd’),title(’Binary Integration Detection Performance’) 
text(l,l,’Pfa’ ,num2str(pfa)) 
legendfOptimal detector’,’binary integration’) 
end 

0. MAIN PROGRAM 

% File name:main.m 
% Title:main program 
% Date of last revision: Mar 10 1995 
% Comments : 

% This m file activates the different programs 

% input varaibles: 

% simpar- File name contains the simulation parameters for the performance comparison 
% Output variables: 

% None 

% Associated Matlab functions: 

% ht.m 

% iht.m 

% alber.m 

% binsnr.m 

% graplot.m 

% Associated Matlab Files: 

% rdata.m 

dat=input(’Enter data file Name 1 ); 

k=menu (Select Required program’,Hough detector’. Hough detector Pd’,Hough detector pfa’,’performance 

comparison’); 

if k==l 

[RT,Q,angnum]=ht(bi,type.dat) 

[rho,tet,mg]=iht(Q,T,dt); 
elseif k==2 
load simpar; 

[pd ]=prdet (pl,p2,sl,s2 ,simlen,bi,type.dat); 
elseif k==3 
load simpar; 

[pfa ]=prfa (si,s2,simlen,bi,type.dat); 
elseif k==4 
load comerpar; 

[sm,pd]=alber(M,N,DC,pfa); 

[SNR]=binpdnr(N,PD,PFA); 
end 
graplot 
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